Merge low coverage visium spots

This document demonstrates how to group and merge visium spots with low coverage using https://github.com/iaaka/visutils The spots are grouped into coarser mesh by seven spots preserving hexagonal layout. All spots from single groups with coverage below certain threshold are merged, their counts are summed and on spot is used as the representative: either the central one (if its coverage below the threshold) or the spot with highest coverage.

First lets install the package, you need to do it only once

devtools::install_github("iaaka/visutils",force = TRUE)
## Using github PAT from envvar GITHUB_PAT
## Downloading GitHub repo iaaka/visutils@HEAD
## 
##   
   checking for file ‘/tmp/RtmpXX0swB/remotes11d8375ad57c/iaaka-visutils-676ab4e/DESCRIPTION’ ...
  
✔  checking for file ‘/tmp/RtmpXX0swB/remotes11d8375ad57c/iaaka-visutils-676ab4e/DESCRIPTION’
## 
  
─  preparing ‘visutils’:
##    checking DESCRIPTION meta-information ...
  
✔  checking DESCRIPTION meta-information
## 
  
─  checking for LF line-endings in source and make files and shell scripts
## 
  
─  checking for empty or unneeded directories
## 
  
─  building ‘visutils_0.0.0.9000.tar.gz’
## 
  
   
## 
## Installing package into '/home/jovyan/R/x86_64-pc-linux-gnu-library/4.0'
## (as 'lib' is unspecified)

Load libraries

library(visutils)
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.0.5
## Attaching SeuratObject
## Attaching sp

Download public dataset from 10x web site

loadVisiumFrom10x("https://cf.10xgenomics.com/samples/spatial-exp/1.3.0","Visium_Mouse_Olfactory_Bulb",outdir="data")
mob = Load10X_Spatial('data/Visium_Mouse_Olfactory_Bulb/')
SpatialFeaturePlot(mob,'nCount_Spatial')

# Quick start everything can be done just in one line

umi.thr = 20000
mobm0 = mergeSpots(mob,getCenters(mob,to.merge = mob$nCount_Spatial<umi.thr))
# lets compare merged object mobm0 with the original one
cex = scaleTo(sqrt(mobm0$nspots),1,sqrt(8),minx = 1,maxx = sqrt(7))*0.9
SpatialFeaturePlot(mob,'nCount_Spatial') + SpatialFeaturePlot(mobm0,'nCount_Spatial',pt.size.factor = cex*2)

Details

Chose coverage threshold

hist(mob$nCount_Spatial,100)

table(mob$nCount_Spatial<20000)
## 
## FALSE  TRUE 
##   457   728

This sample has good coverge and it is not very reasonable to group spots here. However for demonstration purposes I’ll merge all the spots with coverage less than 20000 UMIs. Normally this threshold should be about 500-1000.

groups = getCenters(mob,to.merge = mob$nCount_Spatial<20000)
groups[1:5,]
##                    tissue row col imagerow imagecol  group group.row group.col
## AAACAAGTATCTCCCA-1      1  50 102     2887     6022 49_101        49       101
## AAACCGGGTAGGTACC-1      1  42  28     7078     5238  42_28        42        28
## AAACCGTTCGTCCAGG-1      1  52  42     6284     6223  52_42        52        42
## AAACGAGACGGTTGAT-1      1  35  79     4191     4546  35_79        35        79
## AAACTGCTGGCTCCAA-1      1  45  67     4869     5532  45_67        45        67
mob$group = groups$group
mob@meta.data[1:5,]
##                       orig.ident nCount_Spatial nFeature_Spatial  group
## AAACAAGTATCTCCCA-1 SeuratProject           1055              642 49_101
## AAACCGGGTAGGTACC-1 SeuratProject          31271             6471  42_28
## AAACCGTTCGTCCAGG-1 SeuratProject          31855             6674  52_42
## AAACGAGACGGTTGAT-1 SeuratProject          30391             6803  35_79
## AAACTGCTGGCTCCAA-1 SeuratProject          25989             6610  45_67

Here I color spots that are going to be merged by same color. Note that some of them form hexagons (when all spots have coverage below threshold) and some groups are smaller or even singletons

SpatialPlot(mob,'group',cols = char2col(mob$group,palette = TRUE),pt.size.factor = 2) +NoLegend()

Lets use visutils visualisations since it is a bit cleare (specificall hex one)

par(mfrow=c(1,2),mar=c(1,1,1,1))
plotVisium(mob,mob$group,plot.legend = F)
plotVisium(mob,mob$group,plot.legend = F,type='hex')

Now we have groups, lets use them to merge spots

mobm = mergeSpots(mob,groups)
mobm@meta.data[1:5,]
##                       orig.ident nCount_Spatial nFeature_Spatial group nspots
## GTGGGCTTAGACACAC-1 SeuratProject           1851             1079 28_38      1
## TGTATAACAGATCCTG-1 SeuratProject           8317             3467 28_80      2
## CCTGGCTAGACCCGCC-1 SeuratProject          35378             7453 29_43      5
## CGAGGCTAAATATGGC-1 SeuratProject          31508             7151 29_71      4
## CCACAATGTACGTCTT-1 SeuratProject          21209             5958 29_85      4
##                                                                                                      merged_spots
## GTGGGCTTAGACACAC-1                                                                             GTGGGCTTAGACACAC-1
## TGTATAACAGATCCTG-1                                                          TGCATGTGGTAATCTA-1;TGTATAACAGATCCTG-1
## CCTGGCTAGACCCGCC-1 CCTGGCTAGACCCGCC-1;TCAGGTTCTTTGAGAA-1;TCGAAATTTAGGACCA-1;TGGAGTGATGCGATGA-1;TTAAGCGCCTGACCCA-1
## CGAGGCTAAATATGGC-1                    CGAGAGATGTGAACCT-1;CGAGGCTAAATATGGC-1;GACACAAGGGAAGAAA-1;GGCAATAGTCAATGAG-1
## CCACAATGTACGTCTT-1                    ACCCAACGCCCGTGGC-1;CCACAATGTACGTCTT-1;GAGCTAAGGGCATATC-1;TGAGTGGTCCGTGACG-1

The resultant seurat object contains summed counts, metadata table provides informatin about number of spots merged (nspots) and their identity (merged_spots)

par(mar=c(1,1,1,5))
cex = scaleTo(sqrt(mobm$nspots),1,sqrt(8),minx = 1,maxx = sqrt(7))*0.9
#plotVisium(mobm,mobm$nCount_Spatial,zfun = log1p,cex=cex)
SpatialFeaturePlot(mob,'nCount_Spatial') + SpatialFeaturePlot(mobm,'nCount_Spatial',pt.size.factor = cex*2)